!=======+=========+=========+=========+=========+=========+=========+=$
!       TYPE   : SUBROUTINE
!       PROGRAM: nfourier
!       PURPOSE: fourier-transform the natural-spline interpolation
!                of function Green(tau) 
!                calculate function Green(omega)
!       I/O    :
!       VERSION: 2-16-92
!                29-Nov-95 removal of minimal bug concerning 
!                          DIMENSION of rindata
!       COMMENT: cf J. Stoer R. Bulirsch, Introduction to numerical
!                analysis (Springer, New York, 1980)
!=======+=========+=========+=========+=========+=========+=========+=$
!
       SUBROUTINE nfourier(rindata,coutdata)
       include 'param.dat'
       DIMENSION rindata(L)
       DIMENSION rincopy(L+1),a(L),b(L),c(L),d(L), 
     & u(L+1), q(L+1),xm(L+1)
       COMPLEX*16 coutdata(0:Iwmax),cdummy,explus,ex
       xpi = ACOS(-One)
       delta = Beta/L
       DO i = 1,L
          rincopy(i) = rindata(i)
       ENDDO
       rincopy(L+1) = -1-rindata(1)
       Three = Two+One
       six = Two*Three
     
!c
!c     spline interpolation:  the spline is given by
!c     G(tau) = a(i) + b(i) (tau-tau_i) + c(i) ( )^2 + d(i) ( )^3
!c     The following formulas are taken directly from  Stoer and
!c     Bulirsch p. 102
!c
       q(1) = Zero
       u(1) = Zero
       DO k = 2,L
          p = q(k-1)/Two+Two
          q(k)=-One/Two/p
          u(k)=Three/delta**2*(rincopy(k+1)+rincopy(k-1)-Two*rincopy(k))
          u(k)=(u(k)-u(k-1)/Two)/p
       ENDDO
       XM(L+1) = 0
       DO k = L,1,-1
          XM(k) = q(k)*XM(k+1)+u(k)
       ENDDO
!c
!c     The following formulas are taken directly from  Stoer and
!c     Bulirsch p. 98
!c
       DO j = 1, L
          a(j) = rincopy(j)
          c(j) = XM(j)/Two
          b(j) = (rincopy(j+1)-rincopy(j))/delta - 
     &       (Two*XM(j)+XM(j+1))*delta/6.
          d(j) = (XM(j+1)-XM(j))/(6.*delta)
       ENDDO

!c
!c     The Spline multiplied by the exponential can now be exlicitely
!c     integrated. The following formulas were obtained using
!c     MATHEMATICA
!c
        DO i = 0,Iwmax
           om = (Two*(i)+One)*xpi/Beta
           coutdata(i) = Zero
           DO j = 1,L
              cdummy = ci*om*delta*j
              explus = exp(cdummy)
              cdummy = ci*om*delta*(j-1)
              ex = exp(cdummy)
              coutdata(i) = coutdata(i) + explus*(
     &         ( -six* d(j) )/om**4 + 
     &         ( Two*ci*c(j) + six*delta*ci*d(j)  )/om**3 +
     &         ( b(j)+ Two*delta*c(j)+ three*delta**2*d(j) )/om**2 +
     &         (- ci*a(j) - delta*ci*b(j) - delta**2*ci*c(j) -
     &         delta**3*ci*d(j))/om)
 
              coutdata(i) = coutdata(i) + ex*(
     &        six*d(j)/om**4 - Two*ci*c(j)/om**3 
     &        -b(j)/om**2 + ci*a(j)/om)
           ENDDO
        ENDDO
        END
!
!=======+=========+=========+=========+=========+=========+=========+=$
!       TYPE   : SUBROUTINE
!       PROGRAM: invfourier
!       PURPOSE: inverse fourier transform
!                Greent, Greenw use physical definition
!                Greent(i) = G((i-1)*deltau) for i = 1,...,L
!                Greenw(n) = G(i w_n), for n = 0,L/2-1
!                       w_n = (2*n+1)pi/beta
!                Symmetry property: 
!                G(iw_(-n) = G(iw_(n-1))*
!                coupled to the impurity
!       I/O    :
!       VERSION: 6-16-92
!       COMMENT: 
!=======+=========+=========+=========+=========+=========+=========+=$
!
       SUBROUTINE invfourier(cindata,routdata)
       include 'param.dat'
       DIMENSION routdata(L)
       COMPLEX*16 cindata(0:Iwmax),cdummy
       xpi = ACOS(-One) 
       DO 1 i = 1,L
       routdata(i) = Zero
          tau = (i-1)*beta/real(L)
          DO 2 j = 0,Iwmax
!               om = mod((2*(j)+One)*xpi/Beta*tau,2*xpi)
                om =    ((2*(j)+One)*xpi/Beta*tau)
               cdummy = CMPLX(Zero,One)*om
               dummy = cindata(j)*exp(-cdummy)
         routdata(i) = routdata(i)+Two/beta*dummy
2         CONTINUE
1      CONTINUE
!c
!c     special treatment for tau = 0
!c
       routdata(1) = -One/Two+routdata(1)
       END




!
!=======+=========+=========+=========+=========+=========+=========+=$
! \__ \__  \__ \__ \__ \__ \__ \__ \__ \__ \__ \__ \__ \__ \__ \__ \__ !
! ----------  ----------  ----------  ----------  ----------  ----------
!     TYPE   : SUBROUTINE
!     PROGRAM: FDeSSUN =  = First Derivatives in SU(N) - case
!     PURPOSE: computes first derivatives of GF at G(0) and G(beta) 
!     VERSION: 20-Jun-2003
!     AUTHOR : Viktor Oudovenko
!
! ----------  ----------  ----------  ----------  ----------  ----------
! \__ \__  \__ \__ \__ \__ \__ \__ \__ \__ \__ \__ \__ \__ \__ \__ \__ !
      SUBROUTINE FDeSSUN( omega,eft,G0,G,FM_DOS,UU,occ,docc, FDA1,FDA2)
      include 'param.dat'
	INTEGER io
      REAL*8  DelGF, FM_DOS(nlm,ns),
     &        occ, docc, UU, FDA1(Nd), FDA2(Nd),
     &        omega(nom)
      COMPLEX*16 G(nlm,nom,ns),   G0(nlm,nom,ns)    

      DelGF = 0.0
      DO io = 1,Nom
        DelGF = DelGF+
     &  (ci*omega(io)+eft-1./G0(1,io,1) - FM_DOS(1,1)*1 )*G(1,io,1)
      ENDDO
      DelGF = 2.*DelGF / beta 

! SU(N)  CASE ONLY!
      E = -( eft+0.5*(2*Nlm-1)*UU  - FM_DOS(1,1)*1 )   !???? ef <-> eft
      FD1 = e*(1-occ)-DelGF+(2*Nlm-1)*UU*(occ-docc)
      FD2 = e*occ    +DelGF+(2*Nlm-1)*UU*(    docc)

! manifold FDs
	DO im = 1,Nlm
        DO is = 1,2
	    ims = im + (is-1)*Nlm
          FDA1(ims) = FD1
	    FDA2(ims) = FD2
!          PRINT *,'FDs (SUN)',im,is, E, FDA1(ims), FDA2(ims), DelGF
	  ENDDO
	ENDDO
!
      RETURN
      END
!
!=======+=========+=========+=========+=========+=========+=========+=$
! \__ \__  \__ \__ \__ \__ \__ \__ \__ \__ \__ \__ \__ \__ \__ \__ \__ !
! ----------  ----------  ----------  ----------  ----------  ----------
!     TYPE   : SUBROUTINE
!     PROGRAM: FDeS =  = First Derivatives
!     TYPE   : F77
!     PURPOSE: computes first derivatives of GF  at G(0) and G(beta)
!     VERSION: 20-Jun-2003
!     AUTHOR : Viktor Oudovenko
!
! ----------  ----------  ----------  ----------  ----------  ----------
! \__ \__  \__ \__ \__ \__ \__ \__ \__ \__ \__ \__ \__ \__ \__ \__ \__ !
      SUBROUTINE FDeS( omega, eft,G0,G,FM_DOS,U,zef,nn, FDA1,FDA2 )
      include 'param.dat'
	INTEGER io
      REAL*8  DelGF, FM_DOS(nlm,ns),
     &        zef(nlm,2), U(Nd,Nd), nn(Nd,Nd), FDA1(Nd), FDA2(Nd), 
     &        omega(nom)
      COMPLEX*16 G(nlm,nom,ns),   G0(nlm,nom,ns)    

	DO im = 1,Nlm
        DO ist = 1,2
		is = ist
	    ims = im + (is-1)*Nlm
          IF (Ns .EQ. 1) is = 1

          DelGF = 0.0

          DO io = 1,Nom
            DelGF = DelGF+
     & (ci*omega(io)+eft-1./G0(im,io,is)-FM_DOS(im,is)*1)*G(im,io,is)

          ENDDO
          DelGF = 2.*DelGF / beta 

	    sume = 0.0
	    sumn = 0.0
		sumd = 0.0
 	    DO im1 = 1,Nlm
            DO ist1 = 1,2
		  is1 = ist1
	      ims1 = im1 + (is1-1)*Nlm
            IF (Ns .EQ. 1) is1 = 1

	      IF ( ims1 .EQ. ims ) CYCLE
	      sume = sume + U(ims,ims1)
		  sumn = sumn + U(ims,ims1)*zef(im1,is1)
		  sumd = sumd + U(ims,ims1)*nn(ims,ims1)
	      ENDDO
	    ENDDO

          E = -( eft+0.5*sume - FM_DOS(im,is)*1 )   !??? ef <-> eft
          FDA1(ims) = E*(1.-zef(im,is))-DelGF+(sumn-sumd)
          FDA2(ims) = E*    zef(im,is) +DelGF+(     sumd)

!          PRINT *, 'FDs = ', im,is, E, FDA1(ims), FDA2(ims), DelGF
	  ENDDO
	ENDDO
!
      RETURN
      END
!
!=======+=========+=========+=========+=========+=========+=========+=$
! \__ \__  \__ \__ \__ \__ \__ \__ \__ \__ \__ \__ \__ \__ \__ \__ \__ !
! ----------  ----------  ----------  ----------  ----------  ----------
!     TYPE   : SUBROUTINE
!     PROGRAM: moments
!     TYPE   : F77
!     PURPOSE: computes first and the second moments of Green's Function
!     VERSION: 20-Jun-2003 
!     AUTHOR : Viktor Oudovenko
!
! ----------  ----------  ----------  ----------  ----------  ----------
! \__ \__  \__ \__ \__ \__ \__ \__ \__ \__ \__ \__ \__ \__ \__ \__ \__ !
      SUBROUTINE moments(U, ef, zef, nn, FM_DOS, Vk2, FMA, SMA)
      include 'param.dat'
      REAL*8 ef, E(nlm,2),Vk2(nlm,ns), FM_DOS(nlm,ns),
     &                 zef(nlm,2),U(Nd,Nd), nn(Nd,Nd), FMA(Nd), SMA(Nd)

!	E = -ef
	DO im = 1,Nlm
        DO ist = 1,2
		is = ist
	    ims = im + (is-1)*Nlm
          IF (Ns .EQ. 1) is = 1

 		sum0 = 0.0
	    sum1 = 0.0

		DO im0 = 1,Nlm
		  DO ist0 = 1,2
		    is0 = ist0
 			ims0= im0+(is0-1)*Nlm
              IF (Ns .EQ. 1) is0 = 1
			IF ( ims0 .EQ. ims ) CYCLE
		    sum0 = sum0 + U(ims,ims0)
	        sum1 = sum1 + U(ims,ims0)*zef(im0,is0)
            ENDDO  
	    ENDDO
	    E(im,ist)= -(ef + 0.5*sum0   - FM_DOS(im,is)*1)  ! 
	    FMA(ims) = E(im,ist)+sum1
        ENDDO  
      ENDDO

!
	DO im = 1,Nlm
        DO ist = 1,2
		is = ist
	    ims = im + (is-1)*Nlm
          IF (Ns .EQ. 1) is = 1

	    sum1 = 0.0	
	    sum2 = 0.0
		DO im1 = 1,Nlm
		  DO ist1 = 1,2
		    is1 = ist1
 			ims1= im1+(is1-1)*Nlm
              IF (Ns .EQ. 1) is1 = 1

			IF ( ims1 .NE. ims ) THEN 
			  sum1 = sum1 + 2.*E(im,ist)*U(ims,ims1)*zef(im1,is1)
	        END IF
		
		    DO im2 = 1,Nlm
		      DO ist2 = 1,2
		        is2 = ist2
			    ims2 = im2+(is2-1)*Nlm
                  IF ( Ns .EQ. 1 ) is2 = 1

	            IF ( ims2 .EQ. ims  ) CYCLE
	            IF ( ims2 .EQ. ims1 ) THEN
	             sum2 = sum2 + U(ims,ims1)*U(ims,ims2)*zef(im1,is1)

			    ELSE
				 sum2 = sum2 + U(ims,ims1)*U(ims,ims2)*nn(ims1,ims2)

				ENDIF
	
 		        
                ENDDO  
	        ENDDO
		    	   	        
           ENDDO  
	   ENDDO

	   SMA(ims) = E(im,ist)**2 + sum1 + sum2 + Vk2(im,is)
       ENDDO  
      ENDDO
!
      RETURN
      END
!
!=======+=========+=========+=========+=========+=========+=========+=$
! \__ \__  \__ \__ \__ \__ \__ \__ \__ \__ \__ \__ \__ \__ \__ \__ \__ !
! ----------  ----------  ----------  ----------  ----------  ----------
!     TYPE   : SUBROUTINE
!     PROGRAM: moments
!     TYPE   : F77
!     PURPOSE: computes first and the second moments of Green's Function
!     VERSION: 23-Feb-2001 
!     AUTHOR : Viktor Oudovenko
!
! ----------  ----------  ----------  ----------  ----------  ----------
! \__ \__  \__ \__ \__ \__ \__ \__ \__ \__ \__ \__ \__ \__ \__ \__ \__ !
      SUBROUTINE moment(Nlm, U, ef, occ, docc, FM_DOS, Vk, FM, SM)
!
!    SU(N) version 
      REAL*8  U, ef, occ, docc, Vk, FM, SM,  E, FM_DOS
	INTEGER Nlm, I, NFAC	
!
	IF(Nlm.EQ.1) C2N = 0 
	IF(Nlm.EQ.2) C2N = 3 
	IF(Nlm.EQ.3) C2N = 10 
	IF(Nlm.GE.4) THEN
          C2N = NFAC(2*Nlm-1)/NFAC(2*Nlm-1-2)/2
	ENDIF  !(Nlm.GE.4) 
!	PRINT *, "C2N = ", C2N, NLM, U, Occ, docc, Vk

	E = -( ef+0.5*(2*Nlm-1)*U  - FM_DOS*1)
	FM = e+(2*Nlm-1)*U*occ

	SM = 2.*e*(2*Nlm-1)*U*occ  
	SM = e**2 + SM + 
     &                (2*Nlm-1)*U**2*(occ + 2.*C2N*docc/(2*Nlm-1))+Vk
!
! info:	xmom2 = FM ; xmom3 = -SM	
!
        RETURN
        END
!
! ----------  ----------  ----------  ----------  ----------  ----------
      FUNCTION NFAC(N)
	INTEGER N, NN
	IF(N.GT.16) STOP
	  NN = 1
	DO I = 1, N
		NN = NN*I
	ENDDO

	NFAC = NN

      RETURN
      END
! ----------  ----------  ----------  ----------  ----------  ----------
! \__ \__  \__ \__ \__ \__ \__ \__ \__ \__ \__ \__ \__ \__ \__ \__ \__ !
!       TYPE   : SUBROUTINE
!       PROGRAM: nfourier
!       PURPOSE: fourier-transform the natural-spline interpolation
!                of function Green(tau) calculate function Green(omega)
!       AUTHOR : Viktor Oudovenko
! ----------  ----------  ----------  ----------  ----------  ----------
! \__ \__  \__ \__ \__ \__ \__ \__ \__ \__ \__ \__ \__ \__ \__ \__ \__ !
       SUBROUTINE nfourier1(FD1,FD2,rindata,coutdata)
       include 'param.dat'
       DIMENSION rindata(L)
       DIMENSION rincopy(L+1),a(L),b(L),c(L),d(L), 
     & u(L+1), q(L+1),xm(L+1), W(0:L)
       COMPLEX*16 coutdata(0:Iwmax),cdummy,explus,ex
       xpi = ACOS(-One)
       delta = Beta/L
       DO i = 1,L
          rincopy(i) = rindata(i)
	        W(i-1) = rindata(i)
       ENDDO
       rincopy(L+1) = -1-rindata(1)
	  W(L) = rincopy(L+1)
       three = Two+One
         six = two*three
     
!c
!c     spline interpolation:  the spline is given by
!c     G(tau) = a(i) + b(i) (tau-tau_i) + c(i) ( )^2 + d(i) ( )^3
!c     The following formulas are taken directly from  Stoer and
!c     Bulirsch p. 102
!c
	IF(FD1.GT..99E30) THEN
	 print *,'Natural spline boundary condition is appeared! Check!'
       q(1) = Zero
       u(1) = Zero
	ELSE
	q(1) = -0.5
	u(1) = (3./delta)*( ( rincopy(2)-rincopy(1) )/delta -FD1 )
	ENDIF
       DO k = 2,L
          p = q(k-1)/Two+Two
          q(k)=-One/Two/p
          u(k)=Three/delta**2*(rincopy(k+1)+rincopy(k-1)-Two*rincopy(k))
          u(k)=(u(k)-u(k-1)/Two)/p
       ENDDO
!       XM(L+1) = 0
!
        IF(FD2.GT..99E30) THEN
	 print *,'Natural spline boundary condition is appeared! Check!'
       q(L+1) = Zero
       u(L+1) = Zero
        ELSE
        q(L+1) = 0.5
        u(L+1) = (3./delta)*(FD2 - (rincopy(L+1)-rincopy(L))/delta )
        ENDIF

        XM(L+1) = ( u(L+1)-q(L+1)*u(L))/(q(L+1)*q(L)+1 )  

       DO k = L,1,-1
          XM(k) = q(k)*XM(k+1) + u(k)
       ENDDO
!c
!c     The following formulas are taken directly from  Stoer and
!c     Bulirsch p. 98
!c
       DO j = 1, L
          a(j) = rincopy(j)
          c(j) = XM(j)/Two
          b(j) = (rincopy(j+1)-rincopy(j))/delta - 
     &           (Two*XM(j)+XM(j+1))*delta/6.
          d(j) = (XM(j+1)-XM(j))/(6.*delta)
       ENDDO
!c
!c     The Spline multiplied by the exponential can now be exlicitely
!c     integrated. The following formulas were obtained using
!c     MATHEMATICA
!c
        DO i = 0,Iwmax
           om = (Two*(i)+One)*xpi/Beta
           coutdata(i) = Zero
           DO j = 1,L
              cdummy = ci*om*delta*j
              explus = exp(cdummy)
              cdummy = ci*om*delta*(j-1)
              ex = exp(cdummy)
              coutdata(i) = coutdata(i) + explus*(
     &         ( -six* d(j) )/om**4 + 
     &         ( Two*ci*c(j) + six*delta*ci*d(j)  )/om**3 +
     &         ( b(j)+ Two*delta*c(j)+ three*delta**2*d(j) )/om**2 +
     &         (- ci*a(j) - delta*ci*b(j) - delta**2*ci*c(j) -
     &         delta**3*ci*d(j))/om)
 
              coutdata(i) = coutdata(i) + ex*(
     &        six*d(j)/om**4 - Two*ci*c(j)/om**3 
     &        -b(j)/om**2 + ci*a(j)/om)
           ENDDO
        ENDDO
        END
!
!=======+=========+=========+=========+=========+=========+=========+=$
! \__ \__  \__ \__ \__ \__ \__ \__ \__ \__ \__ \__ \__ \__ \__ \__ \__ !
! ----------  ----------  ----------  ----------  ----------  ----------
!       TYPE   : SUBROUTINE
!       PROGRAM: nfourier2
!       PURPOSE: makes direct fourier-transform of function Green(tau)
!                first interpolating input function and then taken
!                analitical Fourier transformation of obtained function.
!       VERSION: 09-Mar-2000 
!       AUTHOR : VIKTOR OUDOVENKO
!
!  ATT: Provide this information before to  CALL  nfourier2(In,Out,b)
!  1/w    <--> -[G (0)+G (beta)] = -disc
! -1/w^2  <--> -[G'(0)+G'(beta)] = -xmom2
!  1/w^3  <--> -[G"(0)+G"(beta)] = -xmom3
!  Resume:  
!  disc   --> discontinuity: 4 GFs = 1.; 4 Sigma disc = U**2*n(1-n) 
!  & ALWAYS must be > 0 !
!  xmom3  --> sum of the second derivatives : 
!  xmom2  --> sum of the first derivatives  : 
!
!  XMOM2 should have the same     sing as in GF 1/w expantion !
!  XMOM3 should have the opposite sing as in GF 1/w expantion !
!
! ----------  ----------  ----------  ----------  ----------  ----------
! \__ \__  \__ \__ \__ \__ \__ \__ \__ \__ \__ \__ \__ \__ \__ \__ \__ 
      SUBROUTINE nfourier2(disc, xmom2, xmom3, RIndata,COutdata)
      include 'param.dat'
!
       REAL*8 rindata(L), bt,
     &        a(L),b(L),c(L),d(L), xmw(0:L), xm(L+1),
     &        w(0:L),q(0:L),r(0:L),u(0:L),
     &        du(0:L),dn(0:L),p(0:L)
       COMPLEX*16   coutdata(0:Iwmax),cdummy,explus,ex,csum
!
!	PRINT *, 'from Fourier new disc, xmom2,3 = ', disc, xmom2, xmom3
!
       Delta = Beta/L
       DO I = 1,L
          w(I-1) = rindata(I)
       ENDDO
          w(L) = -disc-rindata(1)

! ----------  ----------  ----------  ----------  ----------  ----------
!   Spline interpolation:  
!   G(tau) = a(i) + b(i) (tau-tau_i) + c(i) ( )^2 + d(i) ( )^3
!   For given sum of the first and second derivatives for the END points
!   Written by V. Oudovenko
! ----------  ----------  ----------  ----------  ----------  ----------

      FDsum = xmom2   ! Sum of the first  derivatives at the END points
      SDsum = xmom3   ! Sum of the second derivatives at the END points
!
        yp0 = w(1)-w(0)
        ypn = w(L)-w(L-1)
       ypn2 = w(L)+w(L-2)-2.*w(L-1)
!
       q(0) = 4.D0
       q(1) = 2.D0-0.5D0/q(0)
       r(0) = 0.5D0/q(0)
       p(0) = 2.D0
       p(1) = 2.D0-0.5D0/q(0)
       u(0) = (6.d0/delta)*( (yp0+ypn)/delta -FDsum)+2.*SDsum
      du(0) = u(0)
      dn(1) = 3.D0*ypn2/delta**2-0.5D0*SDsum+0.5D0*u(0)/q(0)
!
      DO k = 1,L-3
           q(k+1) = ( 2.D0-0.25D0/q(k) )
             r(k) = 0.5D0*r(k-1)/q(k)
             u(k) = 3.D0/delta**2*( W(k+1)+W(k-1)-2.*W(k) )
            du(k) = u(k)-0.5D0*du(k-1)/q(k-1)
           p(k+1) = p(k)-(2.d0*r(k)*r(k-1))
          dn(k+1) = dn(k)+(-1)**(k)*du(k)*2.D0*r(k)
      END DO
           u(L-2) = 3.D0/delta**2*( W(L-2+1)+W(L-2-1)-2.*W(L-2) )
          du(L-2) = u(L-2)-0.5D0*du(L-3)/q(L-3)
           r(L-2) = 0.5D0*r(L-3)/q(L-2)
!
         xmw(L-1) = dn(L-2)-du(L-2)*(0.5D0-r(L-3))/q(L-2)
         xmw(L-1) = xmw(L-1)/( p(L-2)-(0.5D0-r(L-3))**2/q(L-2) )
         xmw(L-2) = ( du(L-2)-(0.5D0-r(L-3))*xmw(L-1) )/ q(L-2)
!
      DO k = L-3,1,-1
         XMw(k) = (du(k)+((-1)**k)*r(k-1)*XMw(L-1)-0.5D0*XMw(k+1))/q(k)
      END DO
         XMw(0) = ( u(0)+XMw(L-1)-XMw(1) )/ 4.D0
!
      XMw(L) = SDsum-XMw(0)
      DO k = 1,L+1
          XM(k) = XMw(k-1)
!          PRINT *, k, XM(k)
!
      END DO
! ----------  ----------  ----------  ----------  ----------  ----------
! The following formulas are taken from Stoer and Bulirsch p. 98
! ----------  ----------  ----------  ----------  ----------  ----------

       DO j = 1, L
          a(j) = W(j-1)
          c(j) = XM(j)/2.D0
          b(j) = (W(j)-W(j-1))/delta - (2.D0*XM(j)+XM(j+1))*delta/6.
          d(j) = ( XM(j+1)-XM(j) )/(6.*delta)
       END DO

! ----------  ----------  ----------  ----------  ----------  ----------
!c The spline multiplied by the exponential can be exlicitely integrated. 
!c The following formulas were obtained using MATHEMATICA
! ----------  ----------  ----------  ----------  ----------  ----------

        DO i = 0,Iwmax
           om = (2.D0*(i)+1.D0)*pi/Beta
           coutdata(i) = 0.D0
           DO j = 1,L

              cdummy = ci*om*delta*j
              explus = exp(cdummy)

              cdummy = ci*om*delta*(j-1)
              ex = exp(cdummy)

              coutdata(i) = coutdata(i) + explus*(
     &         ( -6.D0* d(j) )/om**4 + 
     &         ( 2.D0*ci*c(j) + 6.D0*delta*ci*d(j)  )/om**3 +
     &         ( b(j)+ 2.D0*delta*c(j)+ 3.D0*delta**2*d(j) )/om**2 +
     &         (- ci*a(j) - delta*ci*b(j) - delta**2*ci*c(j) -
     &         delta**3*ci*d(j))/om)
!
              coutdata(i) = coutdata(i) + ex*(
     &        6.D0*d(j)/om**4 - 2.D0*ci*c(j)/om**3 
     &        -b(j)/om**2 + ci*a(j)/om)
           ENDDO
        END DO
!
        RETURN
        END

!
! ----------  ----------  ----------  ----------  ----------  ----------
! \__ \__  \__ \__ \__ \__ \__ \__ \__ \__ \__ \__ \__ \__ \__ \__ \__ !
! ----------  ----------  ----------  ----------  ----------  ----------
!       TYPE   : SUBROUTINE
!       PROGRAM: invfourierNEW
!       PURPOSE: INVERSE Fourier Transform :  W -> T
!                Greent, Greenw use physical definition
!                Greent(i) = G((i-1)*deltau) for i = 1,...,L
!                Greenw(n) = G(i w_n), for n = 0,L/2-1   w_n = (2*n+1)pi/beta
!                Symmetry property:  G(iw_(-n) = G(iw_(n-1))
  !     AUTHOR : Viktor Oudovenko
!   tail = 1 with tail (1/omega) proper treatment
!   tail = 0 without  tail substraction
!
! ----------  ----------  ----------  ----------  ----------  ----------
! \__ \__  \__ \__ \__ \__ \__ \__ \__ \__ \__ \__ \__ \__ \__ \__ \__ !
      SUBROUTINE InvFourierNew(cindata,routdata,xmuR)
      include 'param.dat'
!
      REAL*8       routdata(L), xmuR
      COMPLEX*16   cindata(0:Iwmax)
      LOGICAL TAIL	
!
       TAIL = .TRUE.
! 
      IF (TAIL) THEN
        DO i = 0,Iwmax
                   om = (2*(i)+One)*pi/Beta
           cindata(i) = cindata(i)-1.D0/(ci*om-xmuR)
        ENDDO
      ENDIF
!
       DO i = 1,L
          routdata(i) = Zero
                  tau = (i-1)*beta/real(L)
          DO   j = 0,Iwmax
!                       om = mod((2*(j)+One)*pi/Beta*tau,2*pi)
                        om =    ((2*(j)+One)*pi/Beta*tau)
                     dummy = cindata(j)*exp(-ci*om)
               routdata(i) = routdata(i)+Two/beta*dummy
          ENDDO
       ENDDO
!
      IF (TAIL) THEN
       DO i = 1,L
                  tau = (i-1)*beta/real(L)
          routdata(i) = routdata(i) - 
     &        1.*(1. - 1./(1.+exp(xmuR*beta)))*exp(-xmuR*tau)
       ENDDO
      ELSE
!       special treatment for tau = 0
        routdata(1) = -One/Two+routdata(1)
      ENDIF
!
       RETURN
       END

!=======+=========+=========+=========+=========+=========+=========+=$
!
      SUBROUTINE FDfit(Gt,dtau, FD1,FD2)
      include 'param.dat'
      REAL*8  Gt(L), FD1, FD2

              FD1 = (-25.d0*gt(1) + 48.d0*gt(2) - 36.d0*gt(3)+
     &                16.d0*gt(4) -  3.d0*gt(5) ) / 12.d0 / dtau

              FD2 = ( 25.d0*(-1.d0-gt(1)) - 48.d0*gt(L) + 36.d0*gt(L-1)-
     &                16.d0*gt(L-2) + 3.d0*gt(L-3) ) / 12.d0 / dtau      

      RETURN
      END
!
!=======+=========+=========+=========+=========+=========+=========+=$
